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1.  Introduction 

This  semi-aimual  progress  report  contains  a  summary  of  work  accomplished  on  O.  N.  R. 
contract  number  N00014-86-K-0370,  High  Resolution  Radar  Imaging,  during  the  period  from  1 
June  1988  to  30  November  1988. 

p»The  goal  of  this  project  is  to  formulate  and  investigate  new  approaches  for  forming  images 
of  radar  targets  from  spotlight-mode,  delay-doppler  measurements.  These  measurements  could  be 
acquired  with  a  high-resolution  radar-imaging  system  operating  with  an  optical-  or  radio-frequency 
carrier.  Two  approaches  are  under  study.  The  first  is  motivated  by  an  image-reconstruction 
algorithm  used  in  radionuclide  imaging  called  the  confidence-weighted  algorithm',  here,  we  will 
refer  to  this  approach  as  the  chirp-rate  modulation  approach.  The  second  approach  is  based  on 
more  fundamental  principles  which  starts  with  a  mathematical  model  that  accurately  describes  the 
physics  of  an  imaging  radar-system  and  then  uses  statistical-estimation  theory  with  this  model  to 
derive  processing  algorithms;  we  will  refer  to  this  as  the  estimation-theory  approach. 

Work  accomplished  during  the  reporting  period  is  summarized  in  the  following  section. 

2.  Summary  of  Work  Accomplished 

Progress  during  this  reporting  period  has  been  made  on:  a,  extending  the  estimation-theory 
approach  to  include  a  constraint  on  input  signal-to-noise  ratio;  b,  extending  the  estimation -theory 
approach  to  include  a  sieve  constraint  for  stabilizing  image  estimates;  c,  extending  the 
estimation -theory  approach  to  include  a  specular  or  glint  component  in  the  radar-echo  data;  d, 
analyzing  the  performance  of  the  estimation- theory  approach  through  computer  simulations;  and 
c,  modifying  the  chirp-rate  modulation  approach  through  the  introduction  of  the  Wigner-Ville 
distribution.  Some  of  these  areas  are  described  briefly  below  and  more  completely  in  the  appendices. 
A  patent  was  awarded  associated  with  the  chs^-rate  modulation  approach. 

2.1.  Estimation-Theory  Approach  to  Imaging 

During  this  reporting  period,  a  major  effort  has  been  expended  in  implementing  and  con¬ 
ducting  computer  simulations  to  evaluate  the  performance  of  this  imaging  approach  compared  to 
the  conventional  appraoch  based  on  Fourier  transforms.  Preliminary  results  are  reported  in  reference 
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[1],  a  preprint  of  which  is  included  in  Appendix  1.  We  have  extended  the  estimation-theory 
approach  to  include  constraints  on  the  input  signal-to-noise  ratio  and  for  including  sieve  constraints 
to  stabilize  estimated  target  images.  We  have  found  that  such  constraints  can  improve  the  per¬ 
formance  significantly,  as  described  briefly  in  [1]  for  the  SNR  constraint.  A  simulation  program 
is  presently  being  prepared  to  run  equations  of  the  estimation- theory  approach  on  an  Active  Memory 
Technology  Distributed  Array  Processor  having  1024  processors  connected  in  a  mesh  array;  we 
expect  that  this  will  permit  us  to  perform  simulations  with  modest  sized  images  for  performance 
evaluation  studies.  Effort  continued  to  extend  our  model  to  include  specular  components  in  the 
return  signal;  a  brief  status  report  is  contained  in  Appendix  2. 

2.2,  Chirp-Rate  Modulation  Approach  to  Imaging 

A  patent  was  issued  jointly  to  H.  J.  Whitehouse,  of  the  Naval  Ocean  Systems  Center  in  San 
Diego,  and  D.  L.  Snyder  for  the  chirp-rate  modulation  approach  to  imaging  based  on  the  use  of 
the  confidence- weighted  algorithm  [2].  A  copy  of  this  is  in  Appendix  3. 

The  focus  of  our  research  on  the  chirp-rate  modulation  approach  during  the  reporting  period 
has  been  on  modifying  the  image  formation  equations  following  the  introduction  of  the  use  of  the 
Wigner-Ville  distribution  into  the  problem  by  H.  Whitehouse  [3]. 

3.  References 

1.  P.  Moulin,  D.  L.  Snyder,  and  J.  A.  O’Sullivan,  "Maximum-Likelihood  Spectrum  Estimation 
of  Periodic  Processes  from  Noisy  Data,"  submitted  for  presentation  at  the  1989  Conference  on 
Information  Sciences  and  Systems,  Johns  Hopkins  University,  March  1989.  A  preprint  is  in 
App'ifldix  1. 

2.  H.  J.  Whitehouse  and  D.  L.  Snyder,  Imaging  System,  U.S.  Patent  Number  4,768,156,  Aug. 
30,  1988.  A  preprint  is  in  Appendix  2. 

3.  H.  J.  Whitehouse,  "Delay-Doppler  Radar/Sonar  Imaging,”  presented  at  the  Summer  Program 
on  Signal  Processing,  Institute  of  Mathematics  and  Its  Applications,  Univ.  of  Minnesota,  Minneapolis, 
Aug.  1988. 
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ABSTRACT 

We  have  developed  a  new  approach  to  maximum-likelihood  spectrum  estimar 
tion  of  wide-sense  stationary  processes  from  noisy  data.  A  statistical  model  for  the 
data  is  defined:  The  process  whose  spectrum  is  sought  is  wide-sense  stationary, 
periodic  and  Gaussian,  and  its  observations  are  corrupted  by  an  additive  white  noise. 
A  maximum-likelihood  formulation  of  this  problem  has  been  derived,  and  the  equa¬ 
tions  are  solved  numerically  via  the  expectation-maximization  algorithm.  This 
approach  presents  several  attractive  features,  an  important  one  being  that  the  noise 
corrupting  the  observations  is  now  taken  into  account. 

We  present  some  recent  developments  for  this  problem.  The  statistical  perfor¬ 
mance  of  the  new  maximum-likelihood  spectrum  estimator  is  studied  both  theoreti¬ 
cally  and  numerically.  Comparison  with  traditional  estimators  such  as  the  periodo- 
gram  highlight  several  strong  points  of  the  method.  We  also  identify  certain  limita¬ 
tions,  namely  the  instability  of  estimates  for  high  noise  leveb.  These  limitations  can 
be  alleviated  if  a  priori  information  about  the  signal  is  available.  Two  such  problems 
are  discussed  in  which  the  information  at  hand  has  the  form  of  a  constraint  on  the 
input  agnal-to-noise  ratio. 

We  show  how  such  information  can  be  incorporated  in  the  maximum-likelihood 
estimation  procedure.  First  we  assume  the  signal  power  to  be  known.  Theoretical 
issues  of  existence  and  uniqueness  of  the  solution  are  discussed.  We  proceed  with  a 
problem  in  which  the  information  is  leas  complete,  when  only  an  upper-bound  on  the 
signal  power  is  available.  The  statistical  performance  of  both  constrained  estimators 
is  quantitatively  studied. 
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1.  Introduction 

A  promising  approach  to  maximum-likelihood  estimation  of  Toeplitz  constrained  covariamce 
matrices  has  been  proposed  recently  [ij.  Several  further  developments  can  be  considered.  First,  this 
method  also  Z4>plies  to  the  dual  problem  of  spectrum  estimation.  Another  issue  of  interest  is  that  the 
statistical  model  can  account  for  the  presence  of  additive  noise  corrupting  the  observations  and  for 
linear  transformations  of  the  process  whose  covariance  or  spectrum  is  sought.  These  considerations 
have  motivated  a  new  approach  to  high-resolution  delay-doppler  radar  imaging,  where  a  major  goal  is 
to  produce  estimates  of  the  target’s  scattering  function  [2|.  In  the  special  case  of  a  point  target  and  a 
constant  envelope  transmitted  signal,  this  reduces  to  a  spectrum  estimation  problem. 

This  paper  describes  some  recent  developments  for  this  problem.  We  study  the  statistical  per¬ 
formance  of  the  new  maximum-likelihood  spectrum  estimator  both  theoretically  and  numerically. 
Comparison  with  traditional  estimators  such  as  the  periodogram  highlight  several  strong  points  of  the 
method.  We  also  identify  certain  limitations,  namely  the  instability  of  estimates  for  high  noise  levels. 
These  limitations  can  be  alleviated  if  a  priori  information  about  the  signal  is  available.  Two  such 
problems  are  discussed  here  in  which  the  information  at  hand  has  the  form  of  a  constraint  on  the 
input  signal-to-noise  raitio. 

This  p^>er  is  organized  as  follows.  Our  model  is  presented  in  Section  2.  A  maximum-likelibood 
formulation  of  the  problem  is  given  in  Section  3,  and  the  equations  are  solved  via  the  expectation- 
maximization  algorithm.  Section  4  is  devoted  to  a  statistical  performance  analysis  of  this  estimator 
and  a  comparison  with  two  other  methods.  In  Section  5  we  show  how  a  priori  information  on  the  sig¬ 
nal  can  be  incorporated  in  the  maximum-likelihood  estimation  procedure.  First  we  assume  the  signal 
power  to  be  known.  Theoretical  issues  of  existence  and  uniqueness  of  the  solution  are  discussed.  Sec¬ 
tion  5  deals  with  a  le»  complete  knowledge,  where  only  an  upper-bound  on  the  signal  power  is  avail¬ 
able.  The  last  section  is  devoted  to  a  quantitative  study  of  the  statistical  performance  of  both  con¬ 
strued  estimators. 

2.  Model 

The  following  is  derived  from  the  model  presented  in  [l]  for  a  point  target  and  a  constant 
envelope  transmitted  signal.  The  observation  is  an  N-vector  sample  of  a  wide-sense  stationary, 
periodic,  Gaussian  process  corrupted  by  an  additive  noise  : 

r  »  6  +  tn  ,  (2.1) 

where  b  contains  N  consecutive  samples  of  a  zero-mean  periodic  process  with  length  P  >  N,  and  w 
is  an  zero-mean  white  Gaussian  noise  with  variance  No,  uncorrelated  with  b. 

*  This  w<M*k  wis  supported  bjr  contract  number  NOOOli-^K-OSTO  fh)m  the  Office  of  Naval  Research. 


The  perioclicity  aasumption  is  required  to  guarantee  that  the  likelihood  fimction  is  bounded  above; 
therefore,  there  exists  a  maximum-likelihood  estimator  [ij. 

Now  we  define  the  spectral  process  associated  with  to  be  the  DFT  of  one  period  of  6^. 
Assume  that  we  are  interested  in  estimating  only  M  of  the  components  of  this  spectral  P-vector  (1  < 
M  <  P),  the  other  components  being  zero  with  probability  1  ;  let  e  be  this  M-vector.  This  aasumption 
is  introduced  to  deal  with  the  bandlimited  spectra  encountered  in  radar  explications,  which  arise 
because  radar  targets  have  finite  extent  [2].  c  is  a  Gaussian  random  M-vector  with  diagonal  covari¬ 
ance  E,  whose  entries  o^(0>  i  *  0,..,M-1,  are  real  and  positive,  e  and  6  are  related  by  a  linear 
transformation  : 

6  -  r^c  ,  (2.2) 

where  we  have  defined  the  MxN  matrix  P,  consisting  of  the  first  N  rows  and  the  outer  M  columns  of 
the  PxP  DFT  matrix.  The  superscript  f  denotes  the  Hermitian-transpose  operator  on  matrices. 

Our  model  for  the  observations  can  now  be  written  as 

r  -  Plc  -I-  to  .  (2.3) 

The  covariance  matrix  for  r  is  given  by 

Kr  -  Elrrt]  -  PtEP  +  Afo/jv  ,  (2.4) 

where  is  the  AAcAf  identity  matrix. 

3.  Spectrum  Elstimators 

In  this  section  we  introduce  a  maximum-likelihood  spectrum  estimator  for  the  model  (2.3), 
denoted  by  ^ILl.  We  also  define  two  estimators  which  will  be  analyzed  and  compared  to  ours  in  the 
next  section.  The  first  one  is  the  maximum-likelihood  estimator  derived  assuming  noise-free  data, 
denoted  by  MLO  ;  the  second  one  is  the  periodogram. 

3.1.  MLl  Estimator 

From  (2.4),  the  likelihood  function  for  E  is 

I(r,  E)  -  In  det  (P'EP  +  NqIn)  -  54  rf(P^Er  -f-  A'q/n)"'  r  ■  (31) 

Maximizing  the  likelihood  with  respect  to  E  yields  the  necessary  trace  condition  which  the  estimate  E 
must  satisfy  [1,2]: 

Tr  [P(pfEP  +  NolNp(TT^-r^iT  -  No/wXr^SP  -i-  A/o/n)"' P^5E|  -  0  ,  (3.2) 

for  ail  iWkAf  diagonal  matrices  This  trace  condition  is  a  nonlinear  equation  in  E.  Generally  it  can¬ 
not  be  solved  directly  in  closed-form,  so  some  numerical  search  procedure  must  be  implemented.  M 
elegant  solution  is  the  expectation-maximization  (EM)  algorithm  used  in  [1,2|.  An  initial  estimate  E 
is  selected.  At  step  k-(-l  (k  «  0,1,..)  the  estimate  is  updated  according  to 

-  argmax  g(ElE'*’)  (3.3) 

where 

g(El  e‘*’)  «  -54  2*  In  <r*(i)  -  54  "s  .  (3.4) 

i.0  1^  <^(*) 

and 

E(lc(0l®  I  r,E‘‘’]  -  [e'*’  -  E^*'p(PfE^‘'P-l-AfofN)~‘r’E'*’  E‘*’P(pts‘*’p+A/o//v)-‘ 

X  rrt(PfE‘*’P-(-iVo/Ar)'' r^E^*’]  (t, .) .  (3.5) 

This  algorithm  produces  a  sequence  of  estimates 

»*(,)(*■«)  _£(|c(.)|2lr,E‘‘’]  (3.6) 
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having  increasing  likelihood.  It  can  be  shown  that  the  stable  pcwts  this  algorithm  satisfy  the  neces¬ 
sary  trace  condition  for  a  maxinuser  [2].  The  issue  of  uniqueness  is  addressed  in  [3j. 

Special  eaee  :  M—  P  J 

A  closed-form  expression  for  I)  can  be  derived  in  this  special  case: 

-  max(0,r®  -  Ng)  .  (3.7) 

3.2.  MLO  Estimator 

Additive  noise  corrupting  observations  is  usually  not  included  separately  in  approaches  to  spec¬ 
trum  estimation.  This  model  was  assumed  in  [ij.  The  sequence  of  estimates  of  £  is  still  given  by  (3.6) 
and  (3.5),  in  which  we  now  let  Ng  v  0.  We  call  this  the  MLO  estimation.  Clearly  MLO  and  MLl  are 
equivalent  for  noise-free  problems. 

Special  case  :  M 

The  problem  for  which  the  number  of  observations  (N)  is  equal  to  the  number  of  parameters  to 
be  estimated  (M)  is  of  some  practical  interest  It  also  turns  out  that  the  trace  condition  can  be  solved 
in  closed-form  in  this  instance.  The  matrix  F  is  then  invertible,  indicating  the  existence  of  a  one-to 
one  moping  between  r  and  c.  The  MLO  estimator  is  simply 

tf*(0-l(r-frX.]l*,  (3.8) 

where  denotes  (^~‘)^ 

3.3.  Periodogram 

The  periodogram  estimate  of  the  spectrum  is  defined  as  the  (scaled)  magnitude-squared  Fourier 
transform  of  the  N  observations  padded  with  P-N  zeroes  (4).  The  first  M  spectrum  samples  are  then 
given  by 

a^(i)^(P/N)l(rrXi]l‘  (3.9) 

Special  eaee  :  TV—  A/«  P 

When  N  w  M  P,  the  matrix  F  is  equal  to  the  PxP  DFT  matrix  and  the  periodogram  and 
MLO  estimates  are  the  same.  In  this  case,  a  full  period  of  the  process  is  estimated. 

4.  Performance  Analysis 

In  this  section,  we  estimate  £  for  the  model  (2.3)  and  study  the  statistical  performance  of  the 
three  estimators  above.  For  each  method  the  bias  and  variance  are  evaluated,  where 

aas(£)  -  E(£l  -  £  (4.1) 

and 

l^ar(£l-E(£*l-(W-  (“S) 

As  we  shall  see  in  Section  4.3,  the  performance  strongly  depends  upon  the  input  signal  to  noise  ratio 
defined  by 

51VR*.  -  Eo  /  No  ,  (4.3) 

where  Eg  is  the  average  power  of  the  process,  defined  by 

Eo-(l/P)rr(£l.  (4.4) 

From  (4.1)  and  (4.2),  we  derive  the  mean-squared  error  (\SE)  matrix,  defined  by 

MSE  [£]  -  E((£-£)*]  -  Var(E)  (aios[£l)*  .  (4.5) 

The  output  signal  to  noise  ratio  matrix  is  defined  as  follows  : 

5NR^[£1  -  E(£l  (MSE  (£))^  .  (4.6) 
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In  the  following  section,  we  evaluate  the  bias  and  mean-squared  error  for  the  estimators  derived 
in  Section  3.  Whenever  cloeed-form  expressions  for  the  ML  estimates  cannot  be  derived,  computer 
simulations  are  performed.  Typically  3000  realizations  are  generated  for  each  process.  For  a  given  esti¬ 
mator,  (4.1)  and  (4.S)  are  then  estimated  from  the  3000  estimates. 

4.1.  Performance  Analyaia 

Closed  fmn  expressions  for  the  bias  and  mean-squared  error  are  derived  for  MLO  and  the 
periodogram  when  possible.  Simulations  were  carried  out  to  compare  the  performance  of  the  estimar 
tors  for  various  levels  of  input  SNR.  The  performance  was  then  compared  to  the  Cramer-Rao  lower 
bound  for  the  variance  of  unbiased  estimators.  Much  effort  was  made  for  the  special  case  M  =>  N. 
This  provides  insight  into  the  problem  since  the  MLO  equations  can  be  solved  in  closed  form.  The 
choice  of  P  is  free,  so  long  as  P  >  N  [2). 

4.2.  Closed-form  Expressions  for  Estimator  Performance 

(a) MLl 

As  indicated  in  Section  3.1,  no  closed-form  expression  for  the  estimator  is  available,  so  the 
evaluation  of  bias  and  variance  is  obtained  by  computer  simulation. 

(b)  MLO 

Closed-form  expressions  for  MLO  can  be  derived  when  M  =  N.  The  results  are  presented  below. 

Bias 

Combining  (2.3)  and  (3.8),  we  can  write 

a*(,-)=|(c +r-fw)(.)|2.  (4.7) 

Taking  the  expectation  of  (4.7),  we  get 

E(^®(|)1  =  -b  No(rrt)->(,-,i) ,  (4.8) 

which  implies 

Sias[^'(»)!  =  No(rrt)-'(.-,.) .  (4.8) 

The  bias  is  due  to  the  noise  corrupting  the  observations  and  is  proportional  to  its  variance.  The  sensi¬ 
tivity  of  the  bias  to  the  noise  is  determined  by  the  diagonal  entries  of  the  matrix  (rr^)“‘ . 
A'kan-Squared  Error 

Taking  the  expectation  of  (4.7)  squared,  we  obtain 

o  Af-l 

E((^*(.))®1  -  <r«(0  (2-bd,o)  +  No  <^{i)  {  4(rrt)->(,>)  +  2  E  Re[r-HMl"l  ) 

;-o 

N^o  (  2(rrf)-‘(.-,:)2  ^  1^'  r-t(,-,_,32|2  )  .  (4.10a) 

After  some  algebraic  manipulations,  this  expression  can  be  lower-bounded  by 

£^(^*(0)®]  >  2  [«"(.■)  +  No(rrt)->(.-,.T  -  2  ( E[o{i)\  f  .  (4.iob) 

From  (4.8)  and  (4.10),  (4.5)  becomes 

MSE  [^“(01  -  ( £l(^^.-))l  f  +  {No  (rrt)->(M-)  f  ,  (4.11a) 

and 

MSE  [^“(t)i  >  +  2  <72(1)  No  (rrt)-‘(.-.0  -t-  2  (No  (rr’)-^!,*))*  •  (4.11b) 
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(c)  Periodogram 
Bias 

Combining  (2.3)  and  (3.9),  we  write  the  periodogram  estimates  in  the  equivalent  form 

-  (P/N)  I  (rrtc  +  ru;)(t)|2 .  (4.12) 

Taking  the  expectation  of  (4.12),  we  get 

£(^"(t)i  -  (p/A/f  (rrtErrt  +  fVorrtx.-,t) ,  (4.13) 

and 

Bids[ff®(.-)i  -  ( (P/N)  (rrtsrrt)  -  e  x»‘.0  +  (p/^  No  (rrt)(»,i) .  (4.i4) 

The  bias  contains  two  terms.  The  second  is  due  to  the  noise  and  is  proportional  to  Nq.  The  sensitivity 
of  the  bias  to  the  noise  is  deterrained  by  the  diagonal  entries  of  the  matrix  Pr’.  The  other  term  is 
independent  of  Nq.  Even  for  noise-free  observations,  the  periodogram  is  a  biased  estimator  of  E  unless 
rr*  is  the  identity  matrix.  This  would  be  the  case  only  for  N  «  M  —  P  (observation  of  a  ftill  period  of 
the  process)  or  N/M  — *  oo  (infinite  data). 

^/kan-Squared  Error 

Taking  the  expectation  of  (4.12)  squared,  we  obtain 

£t(<rWl.(/«/A/“)[2(j]  ^(j)l(rr»)(».j)l=*P  +  4^o  E  <»"(j)l(rn(«,j)l"(rr’)(M) 

+  2Ag  (rrO(«>f 

+  I  u2(0)  (rrt)(i,  0)*  +  a^(M/2)  (rrt)(i,M/2)2  +  No  s'  r(«,  j)'  I’'  1  •  (4.15a) 

j-n 

This  expression  is  lower-bounded  by 

2  (P/N)2  [  (PP^ErPt  +  NoPPtXi,*)  1®  =  2  (  f  .  (4.15b) 

From  (4.13)  and  (4.15),  (4.5)  becomes 

MSE  [5-®(i))  -  (  El^*(.-))  f  +  [  (P/N)  (PptErPt  -  (N/P)  E  +  NoP^Xi,.)  /  ,  (4.16a) 

and 

ivcE  [a’*(»)]  >  (p'VN®)  ( [(rrfErr^)(»,»f  +  (rp^Errf  -  (n/p)  E)(t,«f  ] 

+  2  No  (PrO(*,»)  (2  PP’EPP^  -  (N/P)  E)(i,.) 

+  2  [  No  (rP’X*'.*)P  )  •  (4.16b) 


4.3.  Simulation  results 

Process  1 

The  first  process  we  consider  is  real  and  has  period  P  =  10.  Its  spectrum  is  symmetric  and  lowpass  (M 
»  5).  Ail  nonzero  spectrum  samples  are  identical ; 

<r2(,)«i,  i=0,..,4. 

The  number  of  observations  is  N  ~  M  ^  5. 

The  noise  variance  Nq  ranges  from  0  to  1.  Figures  1  and  2  show  the  bias  and  SNR,^  for  the  estima¬ 
tors  of  <^(2)  as  a  function  of  SNPi,,  according  to  the  definitions  (4.1),  (4.3),  and  (4.6).  In  the  absence 
of  additive  noise  (  5NR^  — *  oo),  MLl  and  MLO  are  the  same.  Both  are  unbiased  estimators.  The 
periodogram,  however,  is  biased,  and  its  MSE  is  also  larger  than  the  MSE  for  the  ML  estimators. 
When  No  increases  from  0,  the  performance  of  the  estimators  is  roughly  constant  so  tong  as 
remains  above  some  threshold.  For  larger  Nq,  all  three  estimators  exhibit  a  strong  degradation  in 
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perfonnaace.  Comparing  the  thresholds  for  MLO  and  MLl,  we  see  the  tremendous  improvements 
brought  by  taking  the  noise  into  account  in  the  model.  Typically,  for  a  same  SNR,^,  MLl  will  have 
the  same  performance  as  MLO  operating  in  a  20  dB  noisier  environment. 

We  also  notice  that  the  threshold  for  the  periodogram  is  located  at  a  lower  than  for  the 

ML  estimators.  In  Sections  4.2(b)  and  4.2(c),  we  indicated  how  the  sensitivity  of  the  performance  to 
noise  can  be  determined  for  MLO  and  the  periodogram  when  N  »  M.  It  turns  out  that  for  the  flat- 
spectrum  process  considered  here,  the  periodogram  has  a  lower  sensitivity  than  MLO  and  MLl.  This  is 
thought  to  be  due  to  the  smooth  spectrum  used  in  the  simulation. 

Proeeaa  2 

It  has  been  conjectured  that  the  periodogram  does  not  perform  well  for  nonuniform  spectra  [5j. 
This  motivated  our  study  of  a  sharply  peaked  spectrum.  The  process  has  period  P  =  10,  and  a  single 
nonzero  spectrum  component 

<7^(0).  1. 


There  is  just  N  —  M  ~  1  observation. 

Bias  and  SNR^  for  the  estimators  of  0^(0)  are  plotted  as  a  function  of  SlVRi,  in  Figures  3  and 
4.  In  the  atbsence  of  additive  noise,  the  periodogram  is  very  strongly  biased,  and  its  MSE  is  large. 
Furthermore,  in  high-noise  environment  the  periodogram  is  no  longer  more  robust  than  the  ML  esti¬ 
mators.  Clearly,  the  periodogram  is  outperformed  by  MLO  and  MLl.  It  should  also  be  noticed  that  for 
this  process,  the  improvement  of  MLl  over  MLO  is  quite  reduced. 

Computational  Conaiderationa 

The  convergence  rate  of  the  EM  algorithm  depends  on  several  parameters.  The  computation 
time  for  each  iteration  is  of  order  M  The  number  of  iterations  required  for  convergence  of  the 
algorithm  grows  as  M  and  N  increase.  For  VILl,  more  iterations  are  required  as  No  increases,  espe¬ 
cially  in  the  threshold  region  and  beyond.  Typical  figures  are:  for  process  1  with  No  ^  0.1,  30  itera¬ 
tions  are  required  before  the  spectrum  estimates  are  stable;  when  No  ^  1,  300  iterations  must  be  per¬ 
formed.  Our  algorithm  is  implemented  on  a  Maascomp  model  5500.  Running  the  program  on  3000 
realizations  in  the  latter  case  is  typically  completed  in  8  CPU  hours.  We  are  presently  implementing 
these  algorithms  on  a  mesh-connected  1024  processor  (DAP  by  Active  Memory  Technology),  and  we 
expect  a  major  reduction  in  the  time  required  to  produce  estimates. 

4.4.  Crsmer-Rao  Bounds 

In  this  section,  we  study  how  the  MSE  of  the  estimators  considered  so  far  relates  to  the 
Cramer-Rao  bound  on  the  variance.  The  Cramer-Rao  bound  on  the  variance  of  any  unbiased  (UB) 
estimator  of  a  (i)  for  our  model  has  been  found  to  be  [3] 

US-CR[0*(O1  -  (<7*(0  -h  No(rrt)-‘(t,.))2 .  (4.i7) 

From  (4.5)  and  (4.17),  the  MSE  for  an  unbiased  estimator  whose  variance  attains  the  Cramer-Rao 
bound  is  given  by 

AiSE(0^i)l  -  (02(0  iVo(rr0-‘(«,0f  •  (•‘•is) 

-2 

Next  we  state  the  Cramer-Rao  bound  on  the  variance  of  a  biased  (B)  estimator  of  0  (0  : 

B^Rl^i)]  -  C/B-CRl0*(O]  •  (^19) 

From  (4.5),  (4.17)  and  (4.19),  the  MSE  for  a  biased  estimator  reaching  the  Cramer-Rao  bound  is  given 
by 

A4SE(0*(O!  -  (<72(0  +  No(rmi.i)f  (  +  («as[0^Olf  •  (4-20) 

From  the  analytical  expressions  given  for  E(0  (0]  in  Section  4.2,  we  can  now  calculate  the  gradient  of 
E{0*(O]  for  MLO  and  the  periodogram.  Then,  the  minimum  MSE  for  a  biased  estimator  having  the 
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same  bias  as  MLO  and  the  periodogram  is  derived,  and  a  comparison  with  the  actual  MSE  is  made.  No 
ciosedrform  expression  has  been  found  for  MLl. 

MLO 

From  (4.7), 

dn»(.) 

From  (4.9),  (4.20)  and  (4.21),  the  MSE  is  lower-bounded  by 

iW?E„[^*(0]  -  «r*(.-)  -h  2  n*(.)  No  (rrn-‘(t..)  2  [  No  (rrt)-‘(,-,0  ]* 

Periodogram 
From  (4.13), 

-  (p/N)  (rrtrrtx.-,.) . 

Combining  (4.14),  (4.20)  and  (4.23),  the  MSE  is  lower-bounded  by 

M5E„  (^®(0l  -  (P®/N2)  ( ((rrfEiTtXt,*?  +  (rr^Err'  -  (n/p)  sxi.tfi 
+  2  No  (rTtx«,i)  (2  rr^Errt  -  (n/p)  E)(i,.) 

-(-2[No(rrtx.-,.r)- 

Compariaon  o  f  MSE'a  with  Cramer-Rao  bounda 

The  Cramer-Rao  bounds  (4.22)  and  (4.24)  on  the  MSE  are  the  same  as  the  bounds  (4.11b)  and 
(4.16b)  derived  algebraically  from  the  exact  expressions  (4.11a)  and  (4.16a).  Figure  5  shows  how  the 
lower  bounda  compare  with  the  exact  expressions  for  Process  1.  The  actual  MSE’s  are  3-4  dB  above 
their  respective  bounds. 

4.5.  Diseuasion 

The  results  derived  above  suggest  additional  comments  on  a  comparison  between  periodogram 
and  ML  estimators.  Typically  each  component  of  the  gradient  of  E[v  (i)]  given  in  (4.23)  b  much 
smaller  than  unity  (for  the  processes  we  consider),  and  the  Cramer-Rao  bound  on  the  variance  of  the 
periodogram-like  biased  estimator  b  much  smaller  than  the  Cramer-Rao  bound  on  the  variance  of 
unbiased  estimators.  When  the  variance  dominates  the  MSE,  the  periodogram  offers  a  good  MSE  per¬ 
formance.  Thb  was  the  case  for  Process  1.  For  a  less  uniform  spectrum  such  as  the  one  chosen  for 
Process  2,  the  bias  dominates  the  MSE  and  the  periodogram  b  outperformed  by  the  ML  estimators. 

5.  Constrained  maximum-likelihood  estimation 
6.1.  Description  of  the  problem 

An  examination  of  Figures  1-4  suggests  that  MLl  suffers  in  certain  situations.  When  SNRm  b 
low,  the  estimates  are  biased  and  their  variance  b  large.  Although  the  maximum-likelihood  estimator 
b  asymptotically  unbiased  and  efficient,  these  properties  aire  not  guaranteed  in  the  small-sample  prob¬ 
lems  considered  in  Section  4.  Thb  limitation  can  be  alleviated  if  o  priori  knowledge,  such  as  SNE^n,  b 
available.  Since  Ng  b  known,  such  a  constraint  on  the  signal-to-nobe  ratio  can  be  translated  into  a 
constraint  on  the  signal  power  that  must  be  satbfied  by  the  maximum-likelihood  estimates.  Now  we 
show  how  thb  constraint  can  be  incorporated  into  the  HIM  algorithm.  The  constrained  estimates  exbt 
and  are  unique. 

In  Section  5.2,  SNRj,  b  known.  In  Section  5.3,  our  knowledge  b  more  incomplete,  and  only  an 
upper  bound  on  SNRt,  b  available. 


(4.23) 


(4.24) 


(4.21) 

(4.22) 
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5.2.  Known  SiVR^ 

The  eq\iatioDa  for  MLl  presented  in  Section  3.1  can  be  mo^ed  as  follows  to  satisfy  the  con¬ 
straint.  At  each  step  of  the  EM  algorithm,  we  maximize  Q(£  I  E  )  defined  in  (3.4),  subject  to  the 
power  constraint 

&r-i 

E<»^(«)-PE'q-5,  (5.1) 

where  Eq  is  the  signal  power.  The  solution  also  maximizes 

I  tk\  A/— I 

(3(ElE“’)  +  X(E.T2(t)-S).  (5.2) 

where  X  is  a  Lagrange  multiplier.  Taking  the  gradient  of  (5.2)  with  respect  to  E,  we  obtain  a  quadratic 
equation  for  each  spectral  component 

2X<rV«)-<r^(t)  +  C,-0,  (5-3) 


<7.-E(lc(i)l*lr,E‘*’] 

is  calculated  according  to  (3.5).  The  solution  to  (5.3)  is 


;X^0 


X  =  0, 


where  4  is  either  +1  or  -1.  The  equation  for  X  is 


4SX  -  M-  AVl-SCiX  .  (5.5) 

<•0 

In  general  this  nonlinear  equation  in  X  cannot  be  solved  in  closed-form.  Furthermore,  an  ambiguity 
subsists  about  the  choice  of  the  signs  4-  The  latter  problem  is  solved  by  application  of  the  following 
theorem  : 

Theorem 

Assume  that  Co  >  Cj,  i  =*  1,..,M-1.  Then 

(1) 

4--1  ;  j  =  l,..,M-l 

/o-+l  :  5  <  2Co  [M-  S  4VH<^./C’o)' ) 

i-i 

»  —1  ;  else 

(2)  X  is  the  largest  nonzero  solution  of 

(  45X  -  M  +  ^'4VMC,/Co)  )^  =  1  -  8CoX,  for  S  ^  ,  (5.6a) 

«— 1 


M-i 

X  -  0,  for  5-2 


X  is  upper-bounded  by  l/8Co,  and  (5.6a)  can  be  solved  numerically  for  X.  Note  that  the  particular 
case  (5.6b)  is  also  the  solution  to  the  unconstrained  maximization  problem.  Next,  a  (i)(*'“l  is  calcu¬ 
lated  from  (5.4).  The  whole  procedure  is  repeated  at  each  maximization  step  of  the  EM  algorithm. 
Note  that  because  of  the  highly  nonfinear  nature  of  the  problem,  no  analytic  expression  is  available  for 
the  constrained  estimator,  even  in  the  special  case  mentioned  in  (3.7). 


5.3.  Known  upper  bound  on  SNRi, 

In  this  section,  the  o  priori  knowledge  about  has  the  form  of  an  upper  bound.  Our 

approach  parallels  that  of  the  previous  section,  with  the  upper  bound  now  expressed  as  an  inequ^itv 
constraint  on  the  estimated  signal  power.  At  each  step  of  the  EM  algorithm,  we  maximize  Q(£  I  E  ’) 
defined  in  (3.4),  subject  to  the  inequality  constraint 

W— 1 

EAi)<PEo^S,  (5.7) 

iaO 

where  Eq  is  the  upper  bound  on  the  signal  power.  If  the  unconstrained  solution  satisfies  the  upper 
bound,  the  constraint  is  inactive  and  the  estimate  is  given  by  (3.6).  Otherwise,  the  constraint  is 
active,  and  as  in  Section  5.2,  the  solution  is  the  maximizer  of  the  expression  (5.2). 

We  can  expect  the  performance  of  this  estimator  to  be  strongly  conditioned  by  the  choice  of  Eq. 
In  the  limiting  case  £))  — »  oo,  the  constraint  is  always  inactive  and  the  estimator  is  equivalent  to  the 
unconstrained  estimator.  For  the  other  extreme  ease  Eq  — »  0,  the  constraint  is  always  active. 

A.  Simulstion  results 

In  this  section,  we  apply  the  SNR-constrained  estimators  derived  above  to  Process  1,  and  we 
evaluate  numerically  both  their  bias  and  mean-squared  error. 

Figures  6  and  7  give  a  plot  of  the  bias  and  SNR^  for  different  estimators  of  (t^(2)  as  a  function 
of  SNRit,  according  to  the  definitions  (4.1),  (4.3),  and  (4.6).  The  estimators  represented  on  these 
figures  are:  the  two  constrained  estimators  of  Section  5,  respectively  denoted  by  EQ-MLE  and  INEQ- 
^ILE,  and  defined  for  the  (true)  power  constraint  5  ~  5;  the  unconstrained  estimator  MLl  of  Section 
3.1;  and  the  periodogram  PER  of  Section  3.3. 

In  the  absence  of  additive  noise  (  5MZn  MLl  and  EQ-MLE  are  unbiased.  The  periodogram 

and  INEQ-MLE,  however,  are  biased.  For  the  latter,  this  can  be  understood  as  follows.  The  sum  of 
the  M  estimates  is  smaller  or  equal  to  5  a  5,  and  therefore  the  sum  of  all  biases  is  negative.  When  Nq 
increases  from  0,  the  performance  of  the  estimators  is  roughly  constant  so  long  as  5NRn  remains 
above  some  threshold.  For  larger  No,  all  estimators  exhibit  a  degradation  in  performance.  Note  that 
for  the  SNR-constrained  estimators,  each  bias  is  upper-bounded  by  S  —  and  lower-bounded  by 
-  <T®(i).  Comparing  the  performance  in  Figure  2,  we  see  the  favorable  effects  of  incorporating 

SNR  constraints  into  the  problem.  For  low  Nq,  SNR^  is  improved.  This  is  due  to  the  estimates  hav¬ 
ing  a  lower  variance,  which  is  the  dominant  term  in  SNR^.  For  very  noisy  data,  the  performance  of 
the  estimators  is  clearly  improved.  We  can  easily  derive  a  lower  bound  for  SNR,^[<r  (t)]  : 

- ^ - <1. 

max  (5  —  (7®(t ) ,  o®(»)] 

This  bound  is  independent  of  Nq. 

Conelusiona 

In  this  paper,  we  have  described  an  approach  to  spectrum  estimation  from  noisy  data,  based 
upon  a  statistical  model  for  the  observations.  First  we  derive  a  maximum-likelihood  estimator,  and 
evaluate  its  statistical  performance.  A  comparison  is  made  with  two  other  methods  that  do  not  take 
the  additive  noise  into  account.  One  is  the  traditional  periodogram  and  the  other  is  the  maximum- 
likelihood  estimator  derived  for  a  noise-free  model.  It  is  shown  that  in  terms  of  bias  and  MSE,  the 
new  estimator  can  offer  a  better  performance  than  the  latter  ones.  The  improvement  over  the  periodo¬ 
gram  is  noticeable  for  rough  spectra:  The  MSE  was  15  dB  lower  for  the  process  we  considered. 

In  general  however,  the  maximum-likelihood  estimates  are  still  unstable  at  high  noise  leveb.  In 
the  second  step  of  our  study,  we  refine  our  technique  to  improve  the  performance  when  some  side 
information  exists.  We  have  studied  one  such  problem  in  which  some  information  about  the  signal- 
Co-noise  ratio  is  available.  The  performance  for  the  SNR-constrained  estimators  has  been  numerically 
evaluated,  and  compared  with  that  of  the  unconstrained  estimator  and  of  the  periodogram.  The  new 
estimators  perform  significantly  better  than  their  competitors  for  low  SNR*,.  Because  of  the  SNR 
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constraint,  the  estimates  are  not  allowed  to  take  on  the  large  values  that  were  produced  in  the  uncon¬ 
strained  estimation  problem.  This  results  in  the  estimates  having  a  lower  variance.  One  additional 
feature  of  our  approach,  and  an  attractive  one,  is  its  versatility.  Only  a  slight  modification  of  the 
(unconstrained)  algorithm  is  required. 

References 

[1]  M.  L  hfiller,  D.  L.  Snyder,  "The  Role  of  Likelihood  and  Entropy  in  Incomplete-Data  Problems; 
Applications  to  Estimating  Point-Process  Intensities  and  Toeplitz  Constrained  Covariances,” 
F^e.  IEEE,  Vol.75,  No.7,  July  1987. 

[2]  D.  L.  Snyder,  J.  A.  O’Sullivan,  M.  L  ^filler.  The  Use  Of  Maximum-Likelihood  Estimation  For 
Forming  Images  Of  Difiuse  Radar-Targets  From  Delay-Doppler  Data,”  Ih'.h'.k  Trans,  on  Injbr- 
mation  Theory,  to  iq>pear  in  1989. 

[3]  J.  A.  O’Sullivan,  P.  Moulin,  D.  L.  Snyder,  ”Cramep-Rao  Bounds  tor  Constrained  Spectrum  Esti¬ 
mation  with  Application  to  a  Problem  in  Radar  Imaging”,  firoe.  of  the  1988  AUerton  Confer¬ 
ence. 

[4]  W.  Davenport,  W.  Root,  Random  Signals  and  Noise,  McGraw  Hill,  1958. 

[5]  D.  Fuhrmann,  personal  communication. 


-  14  - 


Figure  1.  Biaa{a^{2))  for  Process  1 

MLO  (solid  line),  MLl  (dotted  line), 
periodogram  (dashed  line) 


Figure  3.  Sid«(<7^(0))  for  Process  2 

MLO  (solid  line),  MLl  (dotted  line), 
periodogram  (dashed  line) 


Figure  5.  Comparison  of  SMi,^(a^{2))  and 
Cramer-Rao  bounds  for  Process  1 

MLO  (solid  line),  its  CR  bound  (short- 
dashed  line),  periodogram  (long-dashed 
line),  its  CR  bound  (dotted  line) 


Figure  2.  SMi,^(a^{2))  for  Process  1 

MLO  (solid  line),  MLl  (dotted  line), 
periodogram  (dashed  line) 


Figure  4.  Sl^,^{(r{0))  for  Process  2 

MLO  (solid  line),  MLl  (dotted  line), 
periodogram  (dashed  line) 
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Figure  6.  Biaa(<^(2)) 

MLl  (solid  line),  EQ-MLE  (dotted  line), 
INEQ-MLE  (short-dashed  line),  periodo- 
gram  (long-dashed  line) 


Figure  7.  *  SiVR,^(<T2(2)) 

MLl  (solid  line),  EQ-MLE  (dotted  line), 
INEQ-MLE  (short-d^ed  line),  periodo- 
gram  (long-dashed  line) 

*  Currently  we  are  redrawing  Figure  7  in  which  a 
former  definition  of  the  output  signal-to-noiae  ratio 
waaadopted. 
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K.  E.  Krause 
January  8,  1989 


STATUS  REPORT: 

MAXIMUM  LIKELIHOOD  APPROACH  TO 
SPECULAR  TARGET  IMAGING 


The  following  summarizes  activity  between  June  of  1988  and  January  of  1989 
in  the  statistical  model  formulation  and  imaging  approach  for  the  Maximum 
Likelihood  Estimation  based  imaging  of  delay  and  doppler  spread  specular 
targets. 

At  the  time  of  the  last  report,  two  concepts  were  identified  for  further  analysis 
and  possible  selection  as  the  theoretical  model  to  use  in  the  specular  target 
imaging  problem.  They  were:  (1)  Maximization  of  a  likelihood  function  which 
is  assumed  to  factor  into  a  product  of  identically  structured  likelihood  functions, 
one  for  each  scatterer  in  the  delay-doppler  plane  and  (2)  Application  of  the  EM 
algorithm  to  the  likelihood  function  in  quest  of  a  complete/incomplete  data 
space  formulation  which  would  cause  the  factorization  as  described  in  (1)  above 
to  occur. 

The  consequence  of  the  factorizations  mentioned  is  to  reduce  a  multi¬ 
dimensional  problem  to  the  complexity  of  a  one-dimensional  problem  that  will 
be  solved  many  times,  once  for  each  point  in  the  target  space  grid.  The  idea  in 
investigating  the  EM  algorithm  was  to  provide  a  rigorous  justification  for  this 
likelihood  factorization. 

Recalling  that  each  scatterer  in  the  model  under  investigation  is  assumed  to  be 
characterized  by  a  deterministic  amplitude(to  be  estimated  at  each  point  in  the 
target  space  to  form  the  image)  and  a  random  pha$e(varying  from  known 
exactly  to  uniform  -  and  to  be  integrated  out  in  the  estimation  procedure),  the 
EM  algorithm  was  considered  for  a  complete  data  space  which  consisted  of  a 
scatterer  return  plus  white  noise  constituting  the  signal  for  each  scatterer.  By 
this  construction,  the  signals  for  each  scatterer  were  independent,  hence  the 
likelihood  factorable  in  a  theoretically  rigorous  manner  from  the  start.  The 
remaining  issue  was  to  work  out  the  equations  to  see  what,  if  any, 
computational  complications  might  occur  with  this  formulation.  Limiting  forms 
of  the  random  phase  were  considered  in  some  detail.  Specirically,  phase 
known  exactly  was  first  assumed  and  the  equations  for  the  Expectation  and 
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Maximization  steps  worked  out.  The  result  was  a  one  step  iteration  with 
equation  form  comparable  to  the  solution  that  results  in  the  non-EM  approach. 
Added  in  the  Scheme  of  computational  complexity  was  the  need  to  do  a  linear 
smoothing  in  the  E  step,  which  was  easily  done  in  this  case.  Phase  uniform 
was  considered  next.  The  expectation  to  be  evaluated  took  a  very  complex 
form,  so  it  was  decided  to  look  at  limiting  cases  to  see  how  the  formulation  and 
computational  complexity  would  proceed.  In  considering  a  high  signal  to  noise 
ratio  condition,  the  estimate  could  be  theoretically  calculated  but  would  require 
doing  a  nonlinear  smoothing  to  determine  a  constant  required  in  the  solution. 
The  EM  algorithm  in  this  application  then  seems  to  provide  the  basis  for  a 
rigorously  correct  likelihood  factorization  resulting  in  a  one  step  solution,  but 
with  increased  calculational  complexity. 

In  light  of  this  complexity,  and  with  the  realization  that  any  simulations  and/or 
data  that  the  model  would  be  tested  against  would  likely  use  stepped  frequency 
waveforms(the  complex  envelope  of  which  would  make  factorizations  appear 
reasonable  for  the  integration  times  that  would  be  used),  it  was  decided  to 
proceed  with  a  simulation  to  test  the  first,  non-EM  model  formulation. 
Evaluation  of  results  will  determine  the  necessity  to  proceed  with  the  EM 
approach.  Currently,  the  first  model  concept  is  being  coded  for  imaging  of 
simulated  data  from  simple  generic  targets.  Its  performance  in  comparison 
with  standard  imaging  techniques  will  then  be  studied. 
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[37]  ABSTRACT 

An  improved  imaging  system  has  applications  to  syn¬ 
thetic  aperture  radar,  inverse  synthetic  aperture  radar, 
delay-doppier  radar,  positron-emission  topography  so¬ 
nar.  radionetry  and  other  applications  having  a  target 
image  provided  by  a  series  of  dau  parameterized  by  a 
variable  such  m  9.  A  receiver  structure  includes  a  band- 
paaa  matched-filter,  square  law  envelope-detector,  spe¬ 
cialized  processing  and  convolving  to  produce  the  im¬ 
proved  inuges  irrespective  that  the  radar  signals  have 
practical  side  lobe  structures  and  other  features.  Despite 
the  demands  of  specialiwl  processing  the  architecture 
of  the  algorithm  permits  real-time  implementations. 
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DVUGLNG  SYSTEM 

STATEMENT  OF  GOVERNMENT  INTEREST 

Th*  iavenuon  dcscnbed  hereu  mmy  be  nunuractured  ^ 
and  uacd  by  or  for  die  Govemmeni  of  the  (Jiiitcd  States 
of  America  for  govenunental  purposes  without  the 
paymem  of  any  royalties  thereon  or  therefor. 

BACKGROUND  OF  THE  INVENTION  lo 

This  inveaiioa  rdaies  to  a  method  and  means  for 
improving  an  imaging  system.  In  gteater  particularity  it 
is  ^  a  method  and  means  for  incorporating  the  opera- 
dona  of  resolution  enhancement  and  apriori  tnfonnation 
udliaation  simultaneously  in  the  design  of  an  improved 
imagiag  system.  In  still  greater  particularity  it  is  to 
provide  for  an  improved  method  a^  means  for  improv¬ 
ing  an  imaging  system  that  is  adapted  to  imagiag  in  a 
syathede  aperture  radar,  or  in  an  inverse  synthetic  aper¬ 
ture  radar,  a  radiometer,  a  sonar,  an  electromagnetic  or  ^ 
acoustic  tomographic  system  or  a  related  system  in 
which  there  is  an  interaction  between  the  measurements 
that  are  being  taken  of  physical  phenomenon  and  the 
phenomenon  which  are  being  observed. 

Recently  an  analogy  has  become  recognized  which 
eiists  between  delay-doppier  imaging-radar  systems 
and  tomographic  systems  used  in  clinical  radiology. 
The  analogy  appears  to  hold  the  possibility  of  improv¬ 
ing  radar  inuging  because  the  use  of  matched  filtering 
for  noise  supprestiott  is  suggested  even  by  initial  com-  ^ 
parisons.  an^  more  importandy  because  a  line  of  think¬ 
ing  is  emerging  by  whiidi  new  mathematical  models  for 
the  radar-imaging  problem  might  be  formulated  and 
solved  for  improving  processing.  These  new  models 
account  for  dominant  ^ects  including  noise.  M.  Bern-  35 
feld.  in  his  article  endded  ‘‘Chirp  Doppler  Radar”  Ao- 
atdmtx  /SEE  Kof  7^  No.  4  pp  340-541,  April  19S4. 
made  a  restricted  form  of  this  observation  and  the  re¬ 
stricted  form  also  appears  in  a  different  form  in  the 
work  of  D.  Mensa,  S.  Halevy.  and  G.  Wade  in  their  40 
article  endded  “Coherent  Doppler  Tomography  for 
Microwave  Imagug”  Proettdmgt  IEEE  yoL  71.  No.  2 
pp  234-261.  February  1983.  Both  of  these  articles  draw 
tte  analogy  to  a  tomography  system  wherein  the  data 
available  for  processing  are  in  the  form  of  idealized,  43 
noise-ftee  lin^tegrals  through  the  object  being  im¬ 
aged.  This  type  of  tomography  system  embraces  a  situa¬ 
tion  that  is  well  approzimated  with  X-ray  tomography 
systems  because  X-ray  sources  can  be  highly  collimated 
so  as  to  form  narrow  X-ray  beams  of  high  intensity  that  30 
ate  passed  through  the  object  being  imaged.  Although 
the  analogy  was  articulated  in  these  two  articles,  there 
is  strong  reason  to  believe  that  its  applicability  to  practi¬ 
cal  radar/sonar  signab  of  interest  is  limited  because  the 
ambiguity  functions  normally  associated  with  such  ra-  33 
dar/sonar  signals  do  not  approzhnate  line  distributions 
in  mam  and  thua  do  not  permit  the  evaluation  of  line 
integrals  of  the  scattering  function.  Two  additional 
writinp  dealing  with  frequency-stepped,  chirp-signals 
have  dfscussions  which  ci^fy  this  limitation.  M.  Prick-  60 
ett.  and  C.  Chen  in  “Principles  of  Inverse  Synthetic 
Aperture  Radar  (ISAR)  Imaging.”  IEEE  EASCON 
Eteord,  pp.  340-343.  September  1980  and  M.  Prickett 
and  D.  Wehner  in  “Stepped  Frequency  Target  Imag¬ 
ing”.  Appikattaas  Imag*  Uiidtntanding  and  Spatial  63 
ProeiBing  to  Radar  Signals  for  Automatie  Ship  Classifiea- 
tUm  Workshop.  New  Orleans.  La..  February  1979  dis¬ 
cum  side  lobe  structures  and  other  features  that  cause  a 
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departure  from  idealized  line-integrals  and  the  taut  that 
noise  can  be  non-negiigible  in  some  radar-imaging  situa¬ 
tions.  A  solution  wm  not  evident,  however.  These  arti¬ 
cles  are  included  in  the  Appendix  for  a  reader's  conve¬ 
nience. 

Thus,  there  is  a  continuing  need  in  the  state-of-the-an 
for  a  method  and  means  which  may  permit  the  removal 
of  the  restriction  of  noise-free  line-integrals  so  that  gen¬ 
eral  magnitude  squared  ambiguity  functions  can  be 
accommodated  and  the  recognition  of  the  effects  of 
noise  am  be  developed  for  improved  imaging.  In  this 
discussion  the  ambiguity  function  is  defined  m  the  mag¬ 
nitude  squared  of  the  time-frequency  autocorrelation  or 
crom  correlation  function. 

SUMMARY  OF  THE  INVENTION 

The  present  invention  is  directed  to  providing  a 
means  and  method  for  improving  the  target  imaging 
provided  by  a  series  of  discrete  data  parameterized  by  a 
variable  such  m  an  angle  9  in  an  imagiag  system  which 
receives  data  representative  of  a  physical  phenomena 
and  the  pheonomena  being  observed  and  the  interaction 
of  data  therebetween.  Providing  a  plurality  of  discrete 
data  inputs  each  for  one  of  the  series  of  discrete  data 
enables  a  convolving-processing  in  parallel  to  generate 
two-dimensional  preimage  functions: 

/df,/)-  /  /pdT’/i-dT  -  Tj-fMfdf 

where  fafr'.r)  as  a  function  of  9,r',r  is  the  set  of  avail¬ 
able  data  and  wa  are  chosen  by  the  system  designer. 
Summing  the  two-dimensional  preimage  functions  ena¬ 
bles  a  convolving  the  summed  functions  with  a  circu¬ 
larly  symmetric  function  h  where  h  is  obtained  from  the 
equation: 

where  d  is  the  desired  response  to  a  known  distribution 
p.  In  particular  if  pir'.f^  is  a  two-dimensional  delta 
function,  then  h  is  the  point  spread  function  of  the  imag¬ 
ing  system. 

A  prime  object  of  this  invention  is  to  improve  the 
design  of  an  imaging  system. 

Another  object  is  to  provide  for  an  improved  method 
and  means  for  improving  an  imaging  system  relying  on 
an  interaction  between  the  measurements  taken  of  a 
physical  phenomena  and  the  phenomena  which  is  being 
observed. 

Still  another  object  of  the  invention  is  to  provide  for 
an  improved  imaging  system  relying  upon  new  process¬ 
ing  algorithms  implemented  by  associated  circuitry  that 
provide  improved  visualization  of  targets. 

Another  object  is  to  provide  for  an  improved  method 
and  means  for  imaging  targets  having  specialized  pro¬ 
cessing  for  real  time  implementations. 

Yet  still  another  object  of  the  invention  is  to  provide 
for  an  improved  imaging  system  such  as  in  a  synthetic 
aperture  radar,  an  inverse  synthetic  aperture  radar,  a 
sonar,  an  electrometric  or  acoustic  tomographic  system 
or  related  system  having  a  target  image  provided  by  a 
series  of  data  parameterized  by  a  variable  such  as  angle 
9. 

These  and  other  objects  of  the  invention  will  become 
more  readily  apparent  from  the  ensuing  specification 
and  drawing  when  taken  in  conjunction  with  the  ap¬ 
pended  claims. 
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BRIEF  DESCRIPTION  OF  THE  DRAWINGS 

FIG.  1  is  a  icheniaiic  represencanon  of  a  preferred 
form  of  the  invention. 

FIG.  2  IS  illustrative  for  a  2D  Gaussian  approximaiioa  5 
of  the  2D  ellipsoidal  contour  of  the  ambiguity  function 
of  the  linear  FM  waveform. 

DESCRIPTION  OF  THE  PREFERRED 

EMBODIMENTS  ,0 

First,  a  mathematical  analysis  of  radar  imaging  is 
presented  and  analogized  to  a  related  analysis  concent* 
ing  a  tomographic  imaging  system.  These  two  techno¬ 
logical  discussioas  are  set  forth  to  provide  a  thorough 
appreciation  of  the  salient  features  q(  a  specific  embodi-  IS 
nmnt  of  this  inventive  concept.  It  is  of  course  under¬ 
stood  that  the  analogies  of  these  systems  and  a  conse¬ 
quent  improvement  by  the  inclusion  of  this  invennve 
concept  are  also  appli^lc  to  imaging  synthetic  aper¬ 
ture  radar,  inverse  synthetic  apertme  radar,  radiometric  20 
systems,  sonar  systems  or  othCT  devices  in  which  there 
is  an  interactiott  between  the  measurements  that  are 
being  taken  of  the  physical  phenomenon  and  the  phe¬ 
nomenon  which  are  b^g  observed.  More  specifically 
this  improvement  can  be  incorporated  in  systems  which  25 
enable  the  taking  of  a  series  of  measurements  not  identi¬ 
cally  repeating  the  same  measurements  but  paramater- 
ized  by  a  variable  which  is  generally  designated  as  an 
angle  9.  This  need  not  necessarily  be  a  physical  angle  in 
the  case  of  an  inverse  synthetic  aperture  radar  but  could  30 
be  the  angle  made  by  the  ambiguity  (unction  of  the 
chirp  waveform  relative  to  the  delay  aids  on  each  suc¬ 
cessive  transmission  and  is  the  chirp  rate  of  the  actual 
signal  transmitted.  In  the  case  of  a  radioineter  it  is  the 
real  angle  in  which  the  radiometer  observes  the  scene  of  35 
which  it  is  trying  to  fom  an  image. 

In  other  words,  this  concept  applies  to  a  system  in 
which  there  exists  a  free  parameter  called  an  angle  9 
and  that  a  series  of  measurements  is  made,  each  one  of 
these  measurements  u  a  diflerent  angle  of  Bi, . . .  B*.  It  40 
is  recognized  that  each  measurement  may  itself  be  a 
series  of  submeasurements  in  which  the  angle  is  held 
constant  so  as  to  improve  the  quality  of  the  measure¬ 
ment  at  that  angle  What  is  being  described  isthecombi- 
natkm  of  these  multiple  measurements  and  a  reconstruc-  45 
don  algorithm  which  has  the  capability  of  providing 
high  resolution  and  simultaneously  the  incorporation  of 
aptiori  knowledge  For  example  in  conventional  X-ray 
tomography  an  attenuation  projectiott  can  be  measured 
that  is  sn  integral  of  some  physical  property,  a  scatter-  50 
ing  cron  section  is  measured  in  the  case  of  an  inverse 
synthetic  aperture  radar  and  a  volundnous  flux  projec¬ 
tion  is  measured  in  the  case  of  a  radiometer.  A  set  of 
integrals  is  available  for  these  funedons.  The  fundamen¬ 
tal  theorem  forming  the  subject  matter  of  this  improve-  55 
ment  is  the  radon  inversion  lemma  which  says  that  a 
Fourier  transform  of  this  observation  of  the  one- 
dimensional  projections  is  represented  as  a  slice  through 
the  two-dimensional  Fourier  transform  of  the  distribu¬ 
tion  that  is  trying  to  be  measured.  This,  however,  stand-  60 
ing  by  itself  is  equivalent  to  having  no  additional  infor¬ 
mation  available  and,  therefore,  the  evaluation  of  the 
inverse  radon  transform  tends  to  be  numerically  unsta¬ 
ble  and  is  equivalent  to  doing  numerical  differentiation. 

However,  addition  apriori  information  often  is  avail-  63 
able  For  example  we  know  the  range  resolution  in  the 
case  of  an  inverse  synthetic  aperture  radar,  or  when 
cross-bearing  fixes  are  provided  in  the  case  of  a  radiom- 
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eter.  they  give  an  indication  of  the  approximate  range. 
This  sprion  information  can  be  incorporated  in  pre¬ 
cisely  the  same  way  that  it  is  incorporated  in  the  posi¬ 
tron  emission  tomography  that  is  referenced  herein  and 
disciissrd  in  greater  detail  below.  That  is  to  say.  that 
instead  of  doing  a  back  projection  which  is  a  transfor¬ 
mation  of  a  one-dimensional  data  field  into  a  two-di- 
mensaonal  dam  field  by  spreading  the  one-dimensional 
dam  field  everywhere  pa^el  to  itself  into  two-dimen¬ 
sions,  the  information  is  spread  in  a  region  determined 
by  the  available  apriori  knowledge  This  apriori  knowl¬ 
edge  indicates  a  probability  that  the  information  is  more 
likely  to  be  known  in  one  region  than  it  is  in  another 
region  and  therefore  it  is  not  necessary  to  spread  the 
information  umformly  over  lines  in  the  two-dimensional 
plane  by  back  projectioe  This  operation  is  referred  to 
with  respect  to  position  emission  tomography  as  “confi¬ 
dence  weighting"  but  may  be  interpreted  as  being  a 
weighting  of  the  information  according  to  whatever 
form  of  apriori  knowledge  that  is  available  that  indi¬ 
cates  that  the  information  is  more  likely  to  be  encoun¬ 
tered  in  one  portion  of  the  plane  than  uniformly  along 
lines  (back  projection)  in  the  plane. 

Ra^  imaging  has  been  typified  and  characterized  by 
a  number  of  parameters.  Thm  are  discussed  at  length 
by  H.  L.  Van  Trees  in  his  text  “Dtttetion.  Estimation, 
and  ModuUaion  Thtory:  VoL  3,  John  Wiley  and  Sons, 
New  York,  1971.  The  parameter  p(r,0  is  the  target 
scattering-function  which  is  the  average  reflectivity  as  a 
function  of  delay  r  and  doppler  f,  see  pp.  44g  of  the  Van 
Trees  text  The  parameter  a(r,0  denotes  the  ambiguity 
function  of  the  transmitted  radar-signal,  (page  279  of 
Van  Trees).  In  the  absence  of  noise,  the  output  p(T,f)  of 
a  radar  receiver  consistiag  of  a  bandpass  matched-filter 
(BPMF)  matched  to  the  transmitted  radar  signal)  fol¬ 
lowed  by  a  square-law  envelope-detector  (SLED)  is  the 
convolution  of  the  target  scattered  function  and  the 
ambiguity  function  of  the  transmitted  signal.  Through¬ 
out  this  inventive  concept,  ambiguity  function.  a(r.f) 
refers  to  the  magnitude  squared  ambiguity  function  as 
elucidated  in  somewhat  diflerent  notation  in  the  Van 
Trees  texL  These  expressions  are  set  forth  on  pages  4b2 
and  4d3  of  the  Van  Trees  text  and  form  the  basis  for 

SfiT’/yeT-rz-f^df.  id 

For  the  delay-doppler  radar-imaging  problem  with¬ 
out  noise,  a  sequence  of  target  illuminations  by  chirp- 
FM  signals  is  considered.  Each  of  the  chirp-FM  signals 
has  a  diflerent  chirp  rate.  The  effect  of  changing  the 
chirp  rate  of  a  signal  on  its  ambiguity  ftmetion  is  to 
rotate  the  ambiguity  function  to  an  angle  9  in  the  delay- 
doppler  plane,  (page  291  of  Van  Trees).  This  depen¬ 
dent  is  indicated  in  equation  (I)  by  changing  the  nota¬ 
tion; 


n<ry)-//p<r/)a»<T-r/-/WT'i/<-ii<Ty>  (2) 

where  9  is  determined  by  the  chirp  rate  relative  to  the 
radar  pulse  without  chirp-FM  and  n(r,f)  is  an  undesired, 
naturally  occurring  contaminating  noise  function  which 
is  to  be  minimized  according  to  well  esublished  tech¬ 
niques.  The  noise-free  radar-imaging  problem  arises  in 
the  observation  of  the  output  of  the  BPMF-SLED  re¬ 
ceiver.  p«<r,0  for  a  sequence  of  target  illuminations 
having  diflerent  chirp-FM  rates.  9^9o,  9t,...  9,  and  to 
determine  the  scattering  function  p(r,0. 
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To  eiabontt.  ia  Uie  case  of  (he  chirp  waveform,  (he  Ter-PosossiaB  entiUed  “A  Maihemaiicai  Model  for 
angle  d  not  (he  real  angle  aa  it  would  be  in  the  case  of  Positron  Emisaion  Tomography  Systems  Having  Time- 
trying  to  do  a  triaagulaiion  with  a  radiometer  or  doing  of>F1ight  Measurements'’  /£££  Transaetuna  on  .Vueltar 
(he  actual  physical  measurement  in  (he  positron  emis-  Science  VoL  SS-2t,  pp.  397S-3383.  June  1981.  see  the 
SUM  tomograph.  The  chirp  waveform  angle  is  a  valid  3  appendix. 

parameter,  however.  It  has  been  known  in  the  art  since  The  noise-free  ima^ng  pioblem  of  emission  tomo^a- 
the  fundamental  paper  by  Klauder  in  the  Bell  System  phy  is  to  observe  the  line-of-flight  and  the  time-of-night 
Technical  Journal  in  1980  in  (hat  the  response  of  a  radar  of  the  sequence  of  detected  annihilation  photons,  mod- 
to  a  chirp  waveform  is  parmeterized  by  a  chirp  rate  eled  on  the  average  by  psfr.O  in  equation  (2)  and  to 
number.  This  chirp  rate  number,  the  rate  at  which  the  10  determine  the  two-dimensional  activity  distribution 
frequency  is  changing,  parametetues  the  ambiguity  pfr.IV  Here,  the  parameter  ae(r,0  is  a  known  function 
hinction,  that  is.  its  ability  to  localize  as  a  function  of  determined  by  instrumentation  erron  and  psfr.f)  is  the 
delay  and  doppler.  which  is  inclined  in  (he  delay  dop-  number  of  detected  events  having  a  line-of-flight  with 
pier  plane  at  a  physical  angle  proportioiial  to  the  chi^  angle  9  and  differential  time-of-flight  corresponding  to 
rate,  which  b  the  mathematical  parameter  (hat  de-  11  position  (r,f>  along  the  line-of-flight.  In  a  recent  ezpeh- 
scribes  how  fast  the  chirp  changes.  Stated  in  another  ment  the  dau  has  been  quantized  to  ninety-six  angles 

way,  if  the  chirp  does  not  change  at  all,  (he  angle  is  0  (B/w  |80i/98,  i  wO,  1, . , . ,  9S)  and  to  128-by-l28  posi- 

and  then  one  has  the  ability  to  localize  precisely  in  tions  collecmd  in  an  instrument  being  develop^  at 
doppler  because  there  is  aa  equivaleiicy  to  a  contiauous  Washington  University  and  discussed  by  J.  Blaine.  O. 

sample  of  a  «■«««««<  and  there  is  almost  no  ability  to  20  Ficke,  ti  Hitchens,  and  T.  Holmes  in  their  article  “Data 
resolve  aa  a  ftincdon  of  delay.  As  the  chirp  rate  is  in-  Acquisitioa  Aspects  of  Super-PETT,”  IEEE  Transae- 

creased,  range  and  doppler  are  coupled  together  to  that  thus  on  Svehar  Sehiiet  VoL  iV5-29,  pp.  S4^S47,  Febni- 

a  twodimensional  surface  is  defined  which  relates  the  ary  1982,  see  the  Appendix. 

ability  of  the  waveform  to  resolve  the  target.  By  anal-  The  error  density  aa(r,f)  is  determined  by  both  the 
ogy  to  the  usage  for  optical  imadng  systems,  (he  term  23  physical  size  of  the  crystals  used  in  the  scintillation 
‘'point  spread  function’’  can  be  used  to  describe  this  detectors  (resulting  in  about  a  1-centimeter  uncenainty 
quantity  since  the  ability  of  the  radar  to  resolve  a  single  transverse  of  the  line-of-flight)  and  tne  timing  resolution 
point  target  is  being  dmcribed  and  is  ambiguous  in  the  of  the  electronic  drctiitry  used  to  measure  the  differen- 
sense  that  some  of  the  targets  at  a  near  range  will  be  dal  propogadon-time  (resulting  in  about  a  7s:entimeter 
received  by  the  radar  precisely  the  same  way  for  their  30  spedal  uncertainty  along  the  line-of-flight).  For  present 
doppler  aa  the  targets  at  a  larger  range  will  be  received  systems,  this  density  is  reasonably  modeled  by  a  two-di- 

whh  their  doppler.  That  is,  the  target  point  gives  rise  to  mensional.  eUipdt^y  assymetric  Gaussian-funedon 

a  response  surface  with  a  contour  which  is  an  ellipse  having  its  major  axis  oriented  with  the  line-of-flight  and 

whose  major  axis  is  along  the  line  in  range  delay  dop-  its  minor  axis  oriented  transversely  to  this, 
pier  space  parameterized  by  the  angle  9.  33  For  the  radar-imaging  problem  this  density  corre- 

ReMut  devdopments  ia  positron-emission  tomo-  sponds  to  the  ambiguity  fiiacdon  of  a  radar  pulse  having 
graphic  hnagiag  systems  have  a  relation  analogous  to  aa  envelope  that  is  a  Gaussian  function  and  an  instanta- 

equatkm  (2).  In  these' tomographic  systems  a  positron-  neous  frequoKy  that  is  a  linear  funedon  of  time  and  (he 

emitting  radionueiide  is  introduced  inside  a  padent.  and  phase  which  is  a  quadratic  flmedon  of  time, 

the  resulting  aedvity  is  observed  externally  with  an  <0  From  the  foregoing  for  an  improved  delay-doppler 
array  of  scintillation  detectors  surrounding  the  padent  radar-imaging  system  certain  assumpdons  must  be 
in  a  planar-ring  geometry.  When  a  positron  is  pr^uced  made,  the  firm  of  which  is  that  the  target  is  illuminated 

in  a  radioactive  decay,  it  with  an  electron  by  a  sequence  of  radar  pulses  each  having  a  disdnet 

producing  two  high  energy  photons  that  propagate  in  FM-chirp  rate  correspontog  to  angles  9»9o,  9i . 

oppemu  direcdons  along  a  Uk  In  the  first  system  em-  43  9m  spanning  the  range  from  0*>180’.  A  BPMF-SLED 
ploying  positroa  emission,  the  line-of-flight  of  (he  two  receiver  produces  data  prfr.O  for  9o,  9i.  ■  .  .  9>i  and 
opposheiy  propagating  photons  is  sensed  for  each  de-  quantized  values  of  (r.O-  The  problem  that  remains  is  to 
tected  event.  The  data  attributed  to  these  events  are  estimate  the  target  scattering  function  p(r.f)  using  the 
organized  according  to  their  propagation  angle  and  reiationsbip  stated  in  equation  (2).  For  emission-tomog- 
processed  with  the  same  algorithms  used  in  X-ray  to-  30  raphy  imaging  when  both  dme-of-flight  and  line-of- 
mography.  The  result  is  an  of  the  two-^men-  flight  informadon  are  available,  event  data  is  provided 

sional  spadal  distribution  of  the  radionuclide  within  the  which  is  representadve  of  prfr.O  at  angles  9|, . . 
padent  in  the  plane  of  the  detector  ring.  Recent  devei-  .  ,  9m  spanning  (3-180  and  quantized  to  values  of  (r.O. 
opments  attributed  to  improvements  in  high-speed  elec-  The  measurement-error  density  a«(r,0  is  known.  The 
tronics  and  detector  technology  have  made  it  feasible  to  33  aedvity  distribution  pfr.O  is  to  be  estimated  using  the 
measure  the  useful  accuracy  not  only  of  the  line-of-  reladnship  in  equation  (2).  In  both  cases,  the  delay-dop- 
flight  of  annihilated  photons  but  also  their  difTerendal  pier  radar-imaging  and  the  emission-tomography  imag- 
time-of-flight.  As  a  consequence,  in  the  absence  of  ing,  the  activity  distribution  p<r,0  need  be  estimated 
noise,  the  measurementt  are  in  the  form  of  information  using  (he  relationship  expressed  in  equation  (2). 
of  equation  (2)  with  pfr.f)  beiag  a  two-dimensionai  to  A  number  of  preliminary  considerations  must  be  ex- 
activity  distribution  to  be  iinaged  and  with  aa(r,f)  being  amined  and  defined  to  allow  a  more  thorough  compre- 
the  error  density  associated  with  mesauring  the  location  henskm  of  the  improvement  of  this  inventive  concept, 
of  aa  annihilation  event.  In  this  regard  p(r,f).  aa(r,f).  Pan  of  the  solution  for  improving  the  radar,  tomogra- 
and  pa(r,0  correspond  to  X(x),  p,(x/9).  and  fii[n,P),  phy,  sonar,  or  radiometer  image  lies  in  developing  an 
tcsp^vely,  where  x  and  fs  are  two-dimensional  vec-  63  appropriate  algorithm  for  suiuble  processing  of  an  out- 
tors.  These  parameters  and  their  applicability  to  posi-  put  signal  from  a  known  BPMF-SLED  receiver.  The 
tton  emission  tonwgraphy  systems  are  explained  in  the  algorithm  is  derived  by  applying  sutistical-estimation 
article  by  D.  L.  Snyder,  L.  J.  Thomas,  Jr,  and  M.  M.  theory  to  a  raaihemati^  model  that  accounts  for  (he 
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noise  md  other  efTects  jeen  m  an  emission>(omography 
system  having  lime-of-flight  measuremems.  The  aigo- 
rithm  Tor  soiving  equaiioa  (2)  above  was  proposed  by 
D.  L.  Snyder  et  al  in  their  article  referenced  above.  This 
algonihm  was  evaluated  in  a  later  writing  by  Snyder 
''Some  Noise  Comparnons  of  Data^oUectioa  Arrays 
for  Emission  Tomography-Systems  Having  Time^f- 
Right  Measuremems”  /£££  Tranaaetions  on  Suekar 
Sdinet,  f'oi  iVS-TS,  No.  1.  pp.  1029-1033.  February 
1 982  and  by  Polittc  and  Snyder  in  the  article  ”  A  Simula¬ 
tion  Study  of  Design  Choices  in  the  Implementation  of 
Tune-of-^ght  Recorntruction  Algorithms' 

•agr  IFarfaAgp  on  Tiituhof-flig/u  Tomagnpkji,  Washing¬ 
ton  Univcxsity,  May  198^  published  by  the  igCT-  Com¬ 
puter  Society.  IEEE  Catalog  No.  CHmi-S.  please  see 
these  articles  in  the  Appendix. 

The  noise  was  fou^  to  be  Poisson  distributed  aa 
might  be  expected  because  of  the  quantum  nature  of 
radioactivity  decay,  an  effect  well  modeled  by  a  Pois- 


8 


proach  to  idealize  line-integral  tomography  Examples 
of  weighting  functions  that  might  be  adapted  are: 


itai 


where. 

/Sf<ri/./t/i- 1. 


ri  SSr/1.  /I  Si//la0.  oilwrWHC. 


to 


Here.  ws(r.f)  is  unity  for  delays  and  dopplers  in  a  small 
bin  located  at  r  and  f  in  the  diriay-doppler  plane  and  is 
zero  otherwise,  independea.tly  of  the  sweep  rate  0.  In 


this  ease,  frinf)  equals  pHr.O,  and  the  preimage  is 


IJ  Arji  »  /  ntryws  S  X00ttrj). 

This  choice  of  ws(r.O  might  be  reasonable  if  the  ambi¬ 
guity  fimction  ae(ri)  is  concentrated  about  the  origin 
(rj)m(0ji),  which  requires  a  signal  with  a  large  time- 


concentraiioo  of  the  radkMctive  source.  It  is  argued  in 
the  earlier  referenced  O.  L.  Snyder  et  al  article  that  the 
measured  dam  (that  is,  line-  and  time-of-flight  of  annihi- 
latibn  photons)  are  also  Poisson  distributed,  with  the 
intensity  being  p^rj)  in  equatton  (2).  Mazimum-UIteli- 
hood  esomsoon  is  then  used  to  estiraate  p(rj).  An  ex¬ 
tension  of  this  algorithm  development  is  in  ^ 

later  article  by  O.  L.  Snyder  et  al  entitled  "Image  Re- 
construcdoo  from  List-Mode  Data  in  an  Emisaion  To- 
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mography  System  Having  'nme-of-Right  Measure-  ^  lection  in  tomoeraohy. 

•oubeuM**  twtf  _ area  ^  J 


prciiiiage  is  obtained  simply  by  post  detection  integra¬ 
tion  in  each  dday-dopplm  bin  without  further  process- 
ing. 

Give  w*(rj)  a  value  of  unity  for  values  of  delay  and 
doppler  within  a  narrow  strip  of  width  8  passing 
through  the  origin  of  the  delay-doppler  plane  at  angle  0 
and  W0(r,{)  a  value  of  zero  otherwise  Then  f«(r.O  is  a 
strip  integral,  or  line  integral  for  8  small,  through  the 
data  pt(rj),  which  corresponds  to  unfiltered  baclt-pro- 
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'  /£££  Trmtaetions  on  Sueltar  Scitnct.  VoL  NS- 
20,  No.  3,  pp.  1843-1849,  fune  1983,  see  Appendix.  The 
extended  rigorithm  development  is  said  to  enable  far 
more  aocuraie  reconstructions  at  the  expense  of  greatly 
increased  computation. 

Neglecting  the  eflbets  of  noise  and  statistical  fluctua¬ 
tions  in  the  measurement  data  enable  the  expression  of 
po(rfi  as  the  measurement  described  above.  TIk  im¬ 
proved  imaging  system  10  for  enhanced  radar  imagiag, 
radiometer  imaging,  sonar  imaging,  and  the  Hbr  adapts 
itself  to  the  established  state-of-the-art  and  improves 
thereon,  see  FIG.  L  The  output  pe(r,f)  of  a  BPMF- 
SLED  receiver,  schematically  represented  aa  memory 
10-  is  three-dimensioiial  because  it  is  a  function  of  the 
three  (dependent  variabica  0,  r,  and  f.  The  target  image 
P(r>f)  however,  is  two-dimensionaL  Thus,  a 
three-dimensional  to  two-dimensioiial  transformation  of 
Ps(t4)  is  required  as  part  of  the  improved  processing. 

The  improved  processing  is  accomplisM  in  two 
steps.  The  first  step  is  to  form  a  twodimensional 
"prehnage  array"  20l  This  is  accomplished  by  convolv¬ 
ing  the  dm  psfr.f)  obtained  at  each  FM-chtrp  rate  0 
with  a  weighting  function  wt(r,f)  and  then  summing  the 
results  over  9;  that  is  we  form  the  fiinctions 
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A  confidence  weighting  function  >vs(t'y)««ariT,/)  is 
used  to  fora  the  preimage  as  suggested  from  the  posi- 
tron-emistiott  tomography  experience.  This  corre¬ 
sponds  to  taking  the  value  of  the  BPMF-SLEO  signal 
ps(ri)  shown  in  the  drawings  at  coming  from  measure¬ 
ment  memories  IS  at  each  value  of  delay  and  doppler 
P*i(r4)  •  • .  Psm(t,1)  and  distributing  the  values  over  the 
t^y  doppler  plane  according  to  the  ambiguity  func¬ 
tion  This  approach  is  the  one  now  used  routinely 
in  emission-tomography  systems  having  time-of-flight 
data.  If  the  mathematical  d^elopment  of  the  cited  Sny¬ 
der  et  al  paper  on  mathematical  modeling  carries  over 
the  radar-imaging  problem,  the  choice  of  the  weighting 
function  is  motivated  by  noting  that  the  resulting  frinO 
is  the  maximum-liltelilKxsd  estimate  of  the  delay-dop- 
pier  reflectance  in  the  that  led  to  the  measure¬ 
ment  p0(rj)  assuming  apriori  that  p(r,0  is  uniform. 

The  second  processing  step  in  the  imaging  approach 
relies  on  the  summing  of  the  preimage  outputs  in  a 
summer  2S  and  the  deriving  of  a  target  image  from  the 
preimage  in  a  convolver  30.  Such  a  target  image  resolu¬ 
tion  b  provided  within  a  resolution  function  h(r,0. 
which  defines  a  "desired  image"  according  to 
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/»<r/)»//p»(rV>»(r-r-/-/yr-4r  (3) 

from  which  a  two-ttimensioiial  preimage  ffr.f)  b  de¬ 
rived  according  to 


/ Mr  -  r  /-/Vi<r  /Vr</r 


(4<ll 


-  /  /(KrA«  S  I  Airy). 


(4> 


The  formation  of  thb  preimage  corresponds  to  some 
extent  with  the  back-projection  step  of  the  "unfiltered 
back-projection,  post  two-dimensional  filtering"  sp¬ 


it  has  been  found  that  including  simh  a  resolution 
fiinction  b  important  in  processing  emission-tomogra- 
phy  data  as  a  way  to  trade  off  resolution  and  smoothing 
w  for  nobe  suppression.  A  narrow  two-dimensional,  cir¬ 
cularly  symmetric  Gaussian  resolution-filter  is  used  as 
convolver  30.  Let  d(T,0  denote  the  estimate  of  d(r.n 
obtained  by  processing  the  preimage  fir.f).  Also  let 
0(u.v)  and  F(tt.v)  denote  the  two-dimensional  Fourier 
M  transforms  of  d(T,0  and  fir.O.  respectively.  Thus: 
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•continued 


9 

where  H(u. v)  (see  FIG.  2)  is  the  transform  of  hfr.f)  and 
G<u,v)  is  the  transform  of  the  function  g(r,/)  defined 
according  to 


jtry)  m  I  i/n 


W 

/ 


FIG.  2  is  illustrative  for  the  20  Gaussian  approxtma* 
don  for  the  20  ellipsaidal  contour  of  the  ambiguity  lo 
ftmcdon  of  the  linear  FM  waveform.  This  required 
filter  (unction  is  provided  for  a  goieral  20  contour  in 
accordance  wjth  ettablwhed  techniques. 

Tlw  image  d(ri)  is  obtained  from  6<n.v)  by  a  two-di¬ 
mensional.  inverse  Fourier  transformation.  The  fiinc-  ts 
dons  g(r.i)  and  G<u.v)  are  precomputable  since  they 
depend  o^y  on  the  ambiguity  Amction  and  the 
weighting  function  used  to  form  the  preimage  and  not 
on  the  measured  data. 

For  the  choice  of  the  weighting  (hnctioa:  20 

the  fiinction  g(r.O  is  the  average  over  4  of  the  square  of 
the  ambiguity  fimcdon  aefrj)  is  a  two-dimensional  25 
asymmetric  Gaussiaa  functioa.  and  g(r,0  is  a  Bessel 
funcdoa  The  denvadon  does  not  require  that  Mt(r4)  be 
Gaussiafi.  but  g(r.O  will  usually  need  to  be  evaluated 
numerically  for  practical  ambiguity  Amctiona. 

The  processiag  thusly  described  lends  itself  to  the  30 
radar-imaging  enhancement  and  is  motivated  by  the 
processing  derived  from  a  mathematical  model  for  the 
emission-tomography  imaging  problem.  The  end  result 
is  an  improved  imi^g  for  radar  doppler,  radiometric 
sonar  and  the  like  information  gathe^g  systems.  3] 

The  architecture  suggested  by  the  algorithm  defined 
by  equation  (3)  through  (S)  above  b  similar  to  that  dis¬ 
cussed  in  a  later  article  by  O.  L  Snyder  cntitied  '‘Algo¬ 
rithms  and  Architectures  for  Statuaical  Image  Proc» 
ing  in  Embsioo  Tomography" /lee/ ThneStgiM/Aoecsr-  40 
uif  VIL  VoL  493,  Sod^  of  Photo-Optical  Instrumenta¬ 
tion  Engineen,  pp.  109-111,  19S4,  see  the  Appendix. 
Data  acquired  for  each  doppler  tale  can  be  procested  in 
parallel  and  then  combined  to  form  according  to 
equaiiott  (4)  and  the  processing  in  equation  (3)  required  43 
for  each  doppler  rate  can  be  pipefin^  The  processing 
implemented  in  current  embsioa-tomographs  b  per¬ 
formed  in  the  spatial  rather  than  the  Fourier  domain. 

The  algorithm  has  been  implemented  by  a  computer, 
for  example  a  Perkin-EImm  3242  computer  with  a  30 
floating  point  processor  but  no  array  processor.  Two 
convolutions  and  a  divbiaa  are  required  at  each  stage 
for  each  of  the  data  gatering  angles.  Simultaneous  pro¬ 
cessing  can  be  performed  a^  pipeliaed  for  each  angle 

Obviously  many  modifications  and  variations  of  the  33 
present  invention  are  possible  in  the  light  of  the  above 
teachings  such  as  substituting  magnitude  in  place  of 
magnitude  squared.  It  b  therefore  to  be  understood  that 
within  the  scope  of  the  appended  claims  the  invention 
may  be  practieed  otherudse  than  as  specifically  de-  40 
scribed. 
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Whet  b  claimed  b: 

L  An  apparatus  for  improving  the  target  provided  by 
s  rerws  of  discrete  target  image  data  signab  parameter¬ 
ized  by  a  variable  such  u  an  angle  9  comprising; 
means  for  providing  a  plurality  of  target  image  data 
input  signab  psfrf)  each  for  one  of  the  series  of 
discrete  target  imnge  dam  signals; 
means  coupled  to  a  separate  input  providing  means 
for  processing  the  target  image  dau  input  signals  in 
total  to  generate  separate  two-dimensional  preim- 
age  functions  f(r.f); 

means  for  summing  the  two-dimensional  preimage 
Aumtions; 

means  coupled  to  summing  means  for  convolving  the 
summed  preimage  functions  with  a  response  func¬ 
tion  to  form  an  improved  target  image  signal  ac¬ 
cording  to  the  equation 

^Tj)m  f /Mr-  r/-/Tp<T/)rfrrf/’ 

that  has  equivalence  as 


APPENDIX 
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I.  -Altonibnt  »d  AicMwetunt  lat  SiMiMical  laiwfc 
FrecoMic  w'BaiiauM  m  Touiofnpby- 
by  D.  L.  Sio>dor 


in  Fourier  transform  nouiion  where  H(u.v)  is  the 
Fourier  transform  of  h(r,f);  and 
means  coupled  to  the  convolving  means  for  control¬ 
ling  the  dbpby  of  an  enhanced  image  in  response 
to  Che  improv^  target  image  signal. 
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1  An  apparatus  according  to  cJaiin  1  in  which  each  of 
the  target  inuga  dan  input  signal  providing  means  is 
fabheaiad  to  provide  a  separate  one  of  the  target  image 
dau  input  signab  from  a  separate  one  of  the  series  of 
discrete  target  image  dau  signals  to  be  characteriaed  by 
the  functioa: 

rsiryj-  [  fp(T./>^r-r'/-/XT'/VTV- 
over  0  angicwds . .  ■  d.. 

3.  An  apparatus  according  CO  claim  2  in  which  each  of 
proceming  means  is  fabricated  to  convolve  the  data 
pci(r,f) . . .  with  a  confidence  weighing  Ametion 

wt(Tj)  CO  form  Che  Ametions: 

I 

/dry)"  / 

and  the  summiag  means  sums  over  the  angle  0  to  obtain 
Che  summed  two-dimensional  preimage  Ametion  of: 


ArJi  m  / 
o 

which  optionally  is  expressed  as  the  approxiiiiate  ex- 
pression 

ArJislMrA 

4.  An  apparatus  according  to  claim  3  in  which  the 
convolving  means  is  a  two-dimensional,  circularly  sym- 
metrical  Gaussian  resolution-filter  whose  impulM  re-  ^ 
sponse  is  the  solution  computed  by  a  two-dimensional, 
inverse  Fourier  tnitsformation  expressed  as 


where  a(r.O  is  an  undesired  naturally  occurring 
noise  Ametion  over  0  angtewda . . .  0^ 
processing  in  total  the  plurality  of  series  of  discrete 
target  image  data  input  sign^  to  generate  two-di- 
meniional  preitnage  Ametions  Ur.O; 
summing  the  plurality  of  two-dimensional  preimage 
Ametions;  convolving  the  summed  plurality  of 
twodiatensioiial  preimage  Ametions  into  an  im¬ 
proved  target  image  sig^  corresponding  to  the 
Ametion 


controUittg  the  display  of  the  improved  target  image 
signal  to  provide  an  improved  target  image  on  a 
viewing  screen. 

i.  A  method  according  to  claim  S  in  which  the  step  of 
processing  convoives  each  discrete  target  imagn  dau 
iapttt  sigs^  p$(TJi  with  a  weighting  Aui^on  ws(r,0  to 
be  expressed  as: 


/dryj-  //sslT'/i-dr-  r-/-/Wr 

to  each  obtain  the  two-dimensional  preimage  Ametion 
of: 


for  angle  0^09  ■ . .  0n  which  optionally  is  expressed  as 
the  approxiniate  expression: 


KM  >(!/*)  /  «s<ry)»e(r/WS  AryjaVdrA  ,  -  ...  , 

0  7.  A  method  according  to  claim  6  m  which  the  step  of 

convolving  with  a  circularly  symmetric  response  func- 
that  has  equivalence  as  a  twoKtimensioiial.  circularly  tion  rriies  on  a  two-dimension^  circularly  symmetrical 
symmetrical  Gaussian  resolution-filter  which  satisfies  Gaussian  resolution-filter  computed  by  a  two-dimen- 
tte  equation:  ^  sional  inverse  Fourier  transformation  expressed  as: 


where  G(n,v)  is  the  Fourier  transform  of  g(r,l). 

S.  A  method  of  improving  a  target  image  from  a  43 
series  of  discrete  target  image  data  signals  parameter¬ 
ized  by  a  variable  such  as  an  angle  0  comprising; 
providing  a  plurality  of  discrete  target  image  daa 
input  signals  each  for  one  of  the  series  of  discrete 
target  image  data  signals  and  each  target  image  30 
data  input  signals  expressed  as: 

I  T'/-/yer<r + At./) 


j<ry)  -  (l/v) 


that  has  an  equivalence  as  a  two-dimensional,  circularly 
symmetrical  Gaussian  resolution-filter  which  satisfies 
tte  equation: 

where  G(u,v)  is  the  Fourier  transform  of  g(r,0- 
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